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SUMMARY 


This  report  presents  solutions  of  the  one-  and  two-dimensional 
integral  equations  that  describe  groundwave  propagation.  We  consider 
the  effects  of  (1)  terrain  irregularities  that  are  narrower  than  a 
Fresnel  zone  and  (2)  receiver  elevation.  Our  results  define  the  con¬ 
ditions  under  which  the  simpler  one-dimensional  equation  can  be  used, 
as  well  as  those  that  demand  the  more  complex  two-dimensional  form. 

We  assume  a  frequency  of  100  kHz  throughout,  although  certain  results 
are  easily  scaled  to  other  frequencies. 

It  is  well  known  that  the  one-dimensional  equation  is  invalid 
unless  the  terrain  is  nearly  uniform  across  a  Fresnel  zone.  We 
quantify  the  errors  and  find  that,  for  obstacles  narrower  than  about 
10  km  and  typical  pathlengths,  the  one-dimensional  equation  erron¬ 
eously  predicts  propagation  anomalies  that  (1)  are  independent  of 
width,  and  therefore  too  large,  (2)  diminish  too  slowly  at  long  dis¬ 
tances,  and  (3)  do  not  exhibit  a  diffraction  pattern.  Considerable 
error  can  be  incurred  by  applying  the  one-dimensional  equation  to  moder¬ 
ately  sized  terrain  features.  For  example,  for  a  pathlength  of  500  km, 
that  equation  overstates  by  a  factor  of  4  the  effect  of  an  obstacle 
6  km  in  diameter.  It  cannot  give  accurate  results  unless  the  diameter 
approaches  a  Fresnel  zone  width — which  exceeds  10  km  for  long  propaga¬ 
tion  paths.  Even  for  wide  obstacles,  the  one-dimensional  formulation 
neglects  reflection  and  interference  phenomena. 

The  two-dimensional  integral  equation  has  different  forms  for 
elevated  and  ground-based  receivers.  The  two  forms  become  equivalent 
in  the  limit  of  zero  elevation.  The  stationary-phase  approximation 
can  sometimes  be  used  to  derive  one-dimensional  equations,  which 
are  widely  used  for  most  applications.  This  reduction  in  dimension¬ 
ality  causes  elevation-angle  and  end-point  errors,  both  of  which  are 
associated  with  elevated  receivers. 

Elevation-angle  errors  are  large  unless  the  elevation  is  much 
smaller  than  the  pathlength.  End-point  errors  occur  because  the 
stationary-phase  approximation  improperly  accounts  for  the  region 
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near  the  receiver,  and  are  large  for  elevated  receivers  at  heights 
below  about  one-sixth  of  a  wavelength.  The  one-dimensional  approxi¬ 
mations  for  grounded  and  elevated  receivers  are  therefore  discontinuous 
near  the  ground.  Neither  type  of  error  is  large  for  most  ranges  and 
altitudes  pertaining  to  airborne  receivers.  However,  end-point  errors 
are  significant  below  1  to  2  km,  and  elevation-angle  errors  can  be 
important  for  pathlengths  of  several  tens  of  kilometers  or  less.  In 
practice,  for  receivers  that  are  above  the  ground  but  no  higher  than 
1  or  2  km,  much  of  the  end-point  error  can  be  avoided — without  the 
expense  of  the  two-dimensional  solution — simply  by  using  the  one¬ 
dimensional  solution  for  grounded  receivers. 
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PREFACE 


This  report  continues  Pacific-Sierra  Research  Corporation's 
(PSR)  analysis  of  low-frequency  groundwave  propagation  over  irregular 
terrain.  In  earlier  reports.  Field  and  Allen  derived  validity  cri¬ 
teria  for  application  of  the  one-dimensional  integral  equation  when 
the  terminals  are  on  the  ground;  and  Gayer,  Field,  and  D'Ambrosio 
gave  extensive  numerical  results  based  on  this  equation  for  ground- 
based  terminals  and  generic  terrain  models.^ 

The  present  report  extends  that  earlier  work  by  presenting  and 
comparing  numerical  solutions  of  the  one-  and  two-dimensional  inte¬ 
gral  equations  for  elevated  receivers  and  terrain  features  narrower 
than  a  Fresnel  zone.  In  addition  to  identifying  important  two- 
dimensional  phenomena,  the  results  define  the  conditions  under  which 
the  more  convenient  one-dimensional  equation  can  be  used,  as  well  as 
those  that  demand  the  more  complex  two-dimensional  form. 


E.  Field  and  R.  Allen,  Propagation  of  the  Low-Frequency  Ground- 
wave  over  Nonuniform  Terrain ,  Rome  Air  Development  Center,  RADC-TR- 
78-68,  March  1978. 

+S.  Gayer,  E.  Field,  and  B.  D'Ambrosio,  An  Integral  Equation 
Approach  to  the  Propagation  of  Low-Frequency  Ground  Waves  over  Irreg¬ 
ular  Terrain:  I.  Ground-Based  Terminals ,  Rome  Air  Development  Center, 
RADC-TR-80-334,  November  1980. 
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1 .  INTRODUCTION 


Historically,  two  theoretical  approaches  to  the  analysis  of 
ground-wave  propagation  have  evolved:  the  differential  equation 
method  of  Van  der  Pol,  Norton,  Bremmer,  and  Fock  (Bremmer,  1958]; 
and  the  integral  equation  method  [Hufford,  1952;  Feinberg,  1959]. 

The  differential  equation  yields  virtually  exact  solutions  for  plane 
or  spherical  surfaces  but  is  awkward  for  analyzing  propagation  over 
terrain  whose  conductivity  is  variable  or  whose  topography  is  uneven. 
Under  conditions  of  topographical  irregularity,  numerical  solution 
of  the  integral  equation  is  the  more  practical  approach. 

The  general,  two-dimensional  form  of  the  integral  equation  is 
valid  provided  (1)  the  conductivity  of  the  earth  is  great  enough  to 
justify  applying  impedance  boundary  conditions,  and  (2)  terrain  ir¬ 
regularities  are  not  too  severe.  However,  because  numerical  solu¬ 
tion  of  the  two-dimens ional  integral  equation  is  costly,  its  use  has 
been  limited  to  highly  idealized  irregularities  [de  Jong,  1975],  with 
a  much  simpler,  one-dimensional  approximation  more  commonly  used 
[Johler  and  Berry,  1967;  Johler,  1977;  Gressang  and  Horowitz,  1978]. 

The  one-dimensional  equation  is  derived  by  means  of  a  stationary- 
phase  integration  that  reduces  the  dimensionality  of  the  general  ver¬ 
sion;  being  an  approximation,  it  is  not  valid  for  all  antenna  eleva¬ 
tions  or  terrain  types.  It  is  therefore  necessary  to  define  the 
conditions  under  which  the  more  convenient  one-dimensional  equation 
can  be  used,  as  well  as  those  that  demand  the  more  complex  two- 
dimensional  form. 

Field  and  Allen  [1978]  derived  validity  criteria  for  the  one¬ 
dimensional  equation  as  applied  to  ground-based  terminals;  Gayer, 
Field,  and  D’Ambrosio  [1980]  gave  one-dimensional  numerical  results 
for  ground-based  terminals  and  generic  terrain  models.  The  present 
report  extends  the  earlier  work  by  comparing  numerical  solutions  of 
the  one-  and  two-dimensional  Integral  equations.  For  elevated  re¬ 
ceivers,  we  show  that  the  one-dimensional  results  are  invalid  at 
heights  (1)  Lower  than  a  reduced  wavelength  X/2ir,  or  (2)  so  large 
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as  to  constitute  a  significant  fraction  of  the  pathlength.  For  ground- 
based  receivers,  we  quantify  the  differences  between  the  one-  and  two- 
dimensional  solutions,  finding  serious  disparities  for  terrain  whose 
transverse  irregularities  vary  significantly  across  the  first  Fresnel 
zone.  For  long  paths,  there  are  hence  severe  restrictions  on  the 
terrains  to  which  the  one-dimensional  equation  can  apply. 

Section  II  presents  the  one-  and  two-dimensional  integral  equa¬ 
tions  and  summarizes  the  assumptions  on  which  they  are  based.  Sec¬ 
tion  III  demonstrates  analytically  that  the  one-dimensional  equation 
incorrectly  predicts  discontinuous  fields  near  the  ground,  but  that 
the  two-dimensional  equation  correctly  predicts  continuous  fields; 
it  also  compares  one-  and  two-dimensional  numerical  solutions  as 
functions  of  receiver  elevation  above  a  plane  earth.  Section  IV 
presents  numerical  solutions  for  propagation  over  terrain  that  con¬ 
tains  isolated  inhomogeneities,  and  identifies  two-dimensional  phe¬ 
nomena  not  accounted  for  by  the  one-dimensional  formulation. 

The  combined  results  of  Secs.  Til  and  IV  define  the  regimes  of 
receiver  altitude  and  severity  of  terrain  irregularity  for  which  the 
one-dimensional  treatment  may  be  used.  Section  V  summarizes  the  con¬ 
clusions;  the  appendixes  give  certain  mathematical  details. 
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II.  INTEGRAL  EQUATIONS  FOR  GROUND-BASED  AND  ELEVATED  RECEIVERS 

This  section  summarizes  the  one-  and  two-dimensional  integral 
equations  for  ground-wave  propagation  used  throughout  this  report. 
Each  has  two  forms,  according  to  whether  the  receiver  is  assumed  on 
the  ground  or  elevated.  Detailed  derivations,  having  been  given 
elsewhere  [Hufford,  1952;  Field  and  Allen,  1978],  are  not  repeated. 

We  work  with  the  ground-wave  attenuation  function  W,  defined 
by 


=  2tt|»0  ,  (1) 

where  ip  is  the  vertical  Hertz  potential  and  is  the  Hertz  potential 
in  free  space.  The  function  W  accounts  for  the  earth's  imperfect  con 
ductivity  and  topographic  features.  For  a  flat,  perfectly  conducting 
earth,  W  =  1  and  tp  =  2\p^. 

TWO-DIMENSIONAL  FORMS 

If  surface  impedance  is  small,  terrain  variations  not  too  severe 
and  the  terminals  on  the  ground,  W  satisfies  the  two-dimensional 
integral  equation. 


W(p)  = 


fu2a. 
J  rlr2 


W(Q) 


ik(r ^+r2- 


(2) 


If  the  receiver  is  elevated,  the  following  equation  for  W  applies: 


2W(P)  =  1  + 


ik 

2tt 


0 

*1R2 


W(Q)  e 


ik(r1+R2-RQ 


ng  ^  kR2/  an 


(3) 
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In  Eqs.  (2)  and  (3),  the  transmitter  is  at  the  origin,  and 


k 

n 

n 

g 

P,  P 


da 

Q 


=  free-space  wave  number  w/c, 

=  unit  vector  normal  to  ground, 

=  complex  refractive  index  of  ground, 

=  coordinates  of  ground-based  and  elevated  receivers, 
respectively , 

=  distance  from  transmitter  to  receiver, 

=  element  of  area  on  earth's  surface, 

=  integration  point  on  earth's  surface, 

=  distance  between  transmitter  and  integration  point  Q, 
=  distance  between  Q  and  receiver  points  p,  P. 


Formally,  there  is  no  difference  between  the  definitions  of  p,  r^, 
r^  and  P,  R7,  Rq.  However,  following  Hufford  [1952],  we  use  capital 
Letters  in  Eq .  (3)  to  emphasize  that  the  receiver  is  elevated  above 
the  ground  a  distance  z^. 

Although  Eqs.  (2)  and  (3)  appear  similar,  they  differ  funda¬ 
mentally.  Equation  (2)  is  an  integral  equation  in  that  the  unknown 
function  W(p)  on  the  left  side  also  appears  in  the  integrand  on  the 
right.  Equation  (3)  is  an  integral  formula,  rather  than  equation, 
because  the  unknown  function  W(P)  on  the  left  refers  to  the  attenua¬ 
tion  function  above  the  surface,  whereas  the  integrand  on  the  right 
contains  the  attenuation  function  as  defined  on  the  surface.  Thus, 
to  obtain  W  from  Eq.  (3)  for  an  elevated  receiver,  we  must  first  solve 
F.q .  (2)  to  obtain  W  for  points  on  the  surface,  then  insert  those  values 
into  the  integrand  of  Eq .  (3). 

Equations  (2)  and  (3)  are  the  most  general  forms  of  the  integral 

equation  for  the  ground-wave  attenuation  function  W.  The  quantity 

1/n^  (square  brackets)  represents  the  earth's  surface  impedance  and 

accounts  for  its  imperfect  conductivity.  The  function  3^/90  in 

Eq.  (2)  accounts  for  terrain  shape.  For  terminals  on  a  perfectly 

conducting  flat  earth,  both  i/n  and  9r„/9n  vanish,  and  W  =  1,  as 

R  L 


before. 
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The  function  31^/90  in  Eq.  (3)  accounts  for  both  terrain  shape 
and  receiver  elevation,  and  does  not  vanish  for  a  flat  earth.  A 
careful  limiting  procedure — given  in  Sec.  Ill — is  needed  to  recapture 
the  well-known  result  that  for  a  flat  earth,  W  =  1  if  n  -*■  <»  and 

*0  *  °- 

ONE-DIMENSIONAL  APPROXIMATION 

The  two-dimensional  Eqs.  (2)  and  (3)  are  expensive  to  solve 
numerically.  A  major  simplification  results  if  the  stationary-phase 
method  is  used  for  one  of  the  integrations,  thereby  reducing  Eqs.  (2) 
and  (3)  to  one  dimension.  Details  of  the  reduction  for  ground-based 
receivers  [Eq.  (2)]  are  presented  in  Hufford  [1952]  and  Field  and 
Allen  [1978];  the  one-dimensional  approximation  given  in  those  sources 


W(xQ)  RJ  1 


a("  / 


where  r^,  r^,  and  r^  are  defined  as  above,  the  transmitter  and  re¬ 
ceiver  both  lie  on  the  y  =  0  plane,  and  x  and  x^  are  projections  of 
r^  and  r^  on  the  z  =  0  plane.  Hufford  also  uses  the  stationary-phase 
approximation  to  reduce  Eq.  (3)  (for  elevated  receivers)  to  the  one¬ 
dimensional  form 


2W(xq,  Zq)  —  1  —  e 


— 7Ti/4 


(»)1/2/  dx 


X  W(x 


■,MVVV- 


Both  Eqs.  (4)  and  (5)  are  expressed  in  the  classic  form  derived 
by  Hufford,  which  assumes  that  ^(Rq),  r  ,  and  r2(R9)  nearly  equal 
x^,  x,  and  x^  -  x,  respectively;  that  assumption  requires  terrain  or 
receiver  elevation  to  be  small  compared  with  pathlength.  Field  and 
Allen  [1978J  show  that  somewhat  greater  elevation  can  be  accommodated 
by  substituting — as  we  occasionally  do  below — 


r__J _ i 

Lrir2(rl  +  r2)J 


1/2 


for 


(6) 


in  Kq.  (4)  [and  in  Eq.  (5),  if  we  use  R^,  instead  of  r^,  r2], 
APPLICATION 

Although  some  solutions  have  been  obtained  for  two-dimensional 
equations  equivalent  to  Eqs.  (2)  and  (3)  (see  de  Jong  [1975]  and 
Secs.  Ill  and  TV,  below),  most  calculations  for  irregular  terrain 
iiave  been  based  on  the  far  simpler  Eqs.  (4)  and  (5)  [Johler  and 
Berry,  1967;  Johler,  1977;  Oressang  and  Horowitz,  1978],  It  is 
therefore  important  to  define  the  constraints  on  applicability  im¬ 
posed  by  the  stationary-phase  reduction  of  the  two-dimensional 
equation. 

In  analyzing  the  one-dimensional  equation  for  ground-based  re¬ 
ceivers,  Field  and  Allen  [1978]  show  that  the  reduction  in  dimen¬ 
sionality  is  valid  provided  (1)  the  propagation  path  is  much  longer 
than  a  wavelength,  (2)  terrain  features  have  small  lateral  gradients 
transverse  to  the  direct  propagation  path  (that  is,  gentle  hills  or 
valleys),  and  (3)  terrain  features  vary  only  slightly  across  the 
i irst  Fresnel  zone.  We  show  in  Sec.  IV  that  the  third  condition 
severely  limits  the  applicability  of  the  one-dimensional  equation. 

The  validity  criteria  for  reducing  the  two-dimensional  equation 
tor  elevated  receivers  to  one-dimensional  form  have  not  received 
much  attention  in  the  literature.  We  show  in  the  following  section 
that  Eq .  (5)  is  valid  provided  (1)  the  receiver  is  not  too  close  to 
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the  ground,  and  (2)  the  receiver  height  is  much  less  than  the  path- 
length.  Those  restrictions  are  additional  to  the  ones  already  given 
for  ground-based  receivers. 


T  -  — 
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III.  VALIDITY  REGIME  OF  ONE-DIMENSIONAL  EQUATION 
FOR  ELEVATED  RECEIVERS 


Here,  we  quantify  the  errors  that  result  from  transforming  the 
two-dimensional  equations  to  one-dimensional  form.  We  begin  by 
demonstrating  that  the  two-dimensional  formulation  correctly  repre¬ 
sents  the  continuity  of  the  Hertz  potential  at  low  receiver  eleva¬ 
tions,  whereas  the  approximation  does  not.  Next,  we  show  that  the 
stationary-phase  reduction  from  two  dimensions  to  one  is  inaccurate 
unless  the  elevation  angle  of  the  receiver  is  quite  low.  Finally, 
we  compare  numerical  results  for  both  formulations  and  determine  the 
accuracy  of  the  one-dimensional  equations  as  a  function  of  receiver 
height  and  elevation  angle. 

Because  the  discussion  pertains  solely  to  the  effects  of  the 
receiver's  elevation  above  the  ground,  derivations  and  examples  are 
for  a  smooth-plane  earth.  The  question  of  the  one-dimensional  form's 
ability  to  accurately  represent  propagation  over  irregular  terrain  is 
addressed  in  Sec.  IV. 

CONTINUITY  OF  TWO-DIMENSIONAL  EQUATIONS  NEAR  THE  GROUND 

If  the  receiver  is  on  a  plane  earth,  3^/30  is  identically  zero, 
and  F.q .  (2)  becomes 


W(p)  =  l  + 


ik. 

2tt 


^  elk(rl+r2"r0) 


(7) 


however,  if  the  receiver  is  a  height  z^  above  the  ground, 


SR, 

_ i 

3n 


-z. 


/  2  7  2\l/2  R_  ’ 

\2  +  *o)  2 


(8) 


and  Eq.  (3)  becomes 
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2W(P)  =  1  + 


lk  f  B0  [l  /  i  \  zol  lk(r  «2-k  ) 


Continuity  of  the  fields  requires  that  Eqs.  (7)  and  (9)  become  equal 
in  the  limit  z^  ■*  0.  However,  the  equality  is  by  no  r.eans  immedi¬ 
ately  apparent,  given  the  2  on  the  left  of  Eq.  (9)  ;  ad  the  terms 
proportional  to  Zq/R2  on  tfie 

To  demonstrate  the  low-elevation  equivalence  of  Eqs.  (7)  and 
(9) ,  we  must  recognize  that  even  though 


lim  R„ 


(10a) 


Um  R2  -  r2  , 

V° 


(10b) 


there  is  a  small  region  centered  about  the  subreceiver  point  (near 
r2  =  0)  in  which 


is  undefined  [refer  to  Eq.  (8)],  and  that  a  careful  limiting  procedure 
is  therefore  required.  Subtracting  Eq .  (7)  from  Eq.  (9)  and  applying 
Eqs.  (8)  and  (10),  it  follows  that  continuity  of  the  fields  near  the 
ground  requires 


/'ikro\  f"  f  M<r2'  e)  /  r2z0  \  Ik(rl+r2"r0) 

-""'■hr.  /  d9  /  <*2 -4—  {-JTV * 

'  '  •'O  /o  1  'r2  +  z0' 


2W(P) 


2tt  00 


JO  /*  /  W(r2,  9)  r2zQ  ik(r1+r2-r0) 

2t r  /  d9  /  dr2  /  2  2\3/2  e 

•/0  •/0  1  Ur2  +  zo)  J 


(ID 
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The  r 2  integration  in  Eq.  (14) — found  in  standard  tables — gives  unity, 
and  the  6  integration  obviously  equals  2ir.  Thus,  we  nave  the  correct 
continuity  relation 


-im  1 2W(P)  -  W(p) ]  =  W(p)  .  (15) 

V° 

Two  ooi^ts  should  be  noted  regarding  this  demonstration  of  con¬ 
tinuity  .  First,  the  essential  term  [the  second  on  the  right  in  Eq . 
(11)]  steins  from  (i/kl^)  (BR^/Bn)  in  Eq.  (3).  That  quantity  is  often 
ignored  (e.g.,  Hufford  [1952])  on  the  grounds  that  it  is  a  near-field 
cort..\  tJon  unimportant  for  long  paths — although  ignoring  it  is  clearly 
on  error  for  low-elevation  receivers.  Second,  a  two-dimensional  inte¬ 
gration  is  needed  to  achieve  continuity  of  the  fields. 

DISCONTINUITY  OF  ONE- DIMENSIONAL  EQUATIONS  NEAR  THE  GROUND 

The  discontinuity  of  the  one-dimensional  equations  is  easily 
demonstrated  by  attempting  to  apply  the  limiting  procedure  to  Eqs. 

(4)  and  (5)  rather  than  to  Eqs.  (2)  and  (3).  We  find  that  continuity 
of  the  one-dimensional  equations  requires 


(16) 

in  the  limit  Zq  -*■  0  [just  as  Eq.  (11)  is  required  for  the  two- 
dimensional  case].  Equation  (16)  incorporates  Eq.  (6) — the  substi¬ 
tution  form — and  uses  the  cartesian  notation  P  *  (x^,  z^)  and 
p  *  (Xq,  0),  all  points  lying  in  the  plane  y  =  0. 

It  is  evident  that  Eq.  (16)  cannot  satisfy  the  continuity  re¬ 
quirement  expressed  by  Eq.  (15),  since  the  dependence  on  the  wave 
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number  k  cannot,  be  eliminated  from  either  right-side  term.  For  con¬ 
tinuity  to  apply  at  all  wavelengths,  the  k  dependence  must  disappear 
in  the  limit  z^  -*  0 — as  was  accomplished  in  proceeding  from  Eq .  (11) 
to  Eq.  (15)  for  the  two-dimensional  form. 

The  integrations  in  Eq.  (16)  produce  numerical  inconsistencies 
in  addition  to  the  conceptual  problem  just  mentioned.  To  illustrate, 
we  proceed  as  before  [Eqs.  (11)  to  (15)],  arguing  that  l/R^  peaks 
sharply  about  x  =  x^  and  factoring  all  terms  outside  the  integrals. 
After  some  rearrangement  and  a  change  of  variable,  we  have  the  fol¬ 
lowing  one-dimensional  relation  to  replace  Eq.  (14): 


The  integrals  can  be  evaluated  in  terms  of  the  incomplete  beta  func¬ 
tion  [Gradshteyn  and  Kyzhik,  1965].  Here,  we  simply  note  that  the 

1/2 

first  integral  is  proportional  to  z  and  vanishes  as  z  -*•  0,  and 

1/2  <; 

that  the  second  is  proportional  to  1/z^  and  diverges  as  z^  -*•  0. 

We  also  recall  that  the  second  integral  results  from  the  (often  ig¬ 
nored)  near-field  term  in  Eq.  (5) — ( i/kR2) (3R2/3n) . 

If  near-fields  are  retained,  the  one-dimensional  equation  for 
elevated  receivers  thus  gives  an  undefined  result  as  altitude  ap¬ 
proaches  zero.  If  near-fields  are  ignored,  however,  the  right  side 
of  Eq.  (17)  vanishes  and  continuity,  as  expressed  by  Eq.  (15),  is 
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violated  by  a  factor  of  2.  In  either  case,  the  one-dimensional  equa¬ 
tion  for  elevated  receivers  is  clearly  in  error  if  the  altitude  is 
less  than  what  we  define  below. 

END-POINT  AND  ELEVATION-ANGLE  ERRORS  IN  ONE-DIMENSIONAL  EQUATIONS 


The  stationary-phase  reduction  of  the  two-dimensional  equations 
is  summarized  in  Appendix  A.  For  receivers  on  a  plane  earth,  the 
reduction  in  dimensionality  is  accurate  provided  kxR  »  1.  For  ele¬ 
vated  receivers,  the  reduction  is  subject  to  two  additional  constraints, 
which  we  call  "end-point"  and  "elevation-angle"  errors. 

The  end-point  error  arises  because  (as  discussed  above)  if  re¬ 
ceiver  elevation  is  low,  the  term  (i/kR,,)  (3R2/9n)  in  Eq.  (3)  peaks 
very  sharply  near  the  end  of  the  propagation  path.  The  end-point 
term,  which  is  not  accounted  for  in  the  standard  stationary-phase 
treatment,  must  undergo  a  full  two-dimensional  integration  over  a 
small  region  centered  about  the  subreceiver  point.  The  correction — 

details  of  which  are  given  in  Appendix  A — that  must  be  added  to  the 

•j* 

right  side  of  the  one-dimensional  Eq.  (5)  is  as  follows: 


end-point  correction 


W(xR,  0)  , 


W(xQ,  0)/kzQ  , 


if  2ttzr/A  <<  1  ; 


if  2ttzr/A  »  1  . 


Equation  (18)  shows  that  the  end-point  correction  will  cancel 
half  the  left  side  of  Eq .  (5)  if  2fzr/A  «  1,  thereby  removing  the 
factor  of  2  and  causing  Eqs.  (4)  and  (5)  to  yield  the  same  value  for 
W  near  the  ground.  Omitting  the  correction  under  the  same  conditions 
will  cause  significant  errors.  However,  Eq.  (19)  shows  that  as 


Although  corrections  to  these  errors  are  derived  for  a  plane 
earth,  they  pertain  solely  to  receiver  elevation  and  should  apply 
even  for  irregular  terrain. 

To  avoid  counting  the  end-point  region  twice,  and  to  avoid  the 
low-elevation  divergence  discussed  in  connection  with  Eq.  (17),  we 
must  also  exclude  the  region  surrounding  x  Xq  when  integrating  the 
terms  in  Eq.  (5)  that  contain  SR^/Sn. 
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receiver  elevation  increases,  the  end-point  error  in  the  one-dimen¬ 
sional  equation  becomes  smaller — being  in  fact  negligible  if 
2TT20/X  »  1. 

Certain  of  the  approximations  used  to  transform  Eq .  (3)  into 
Eq.  (5)  break  down  (see  Appendix  A)  if  the  receiver  elevation  angle 
z^/Xq  becomes  too  large.  That  constraint  on  elevation  angle  can  be 
combined  with  Eqs.  (18)  and  (19)  to  give  the  following  criteria  for 
the  domain  of  validity  of  the  one-dimensional  equation: 


_X_ 

2u 


« 


(20) 


The  lower  limit  on  receiver  height  is  the  end-point  error;  the  upper 
limit  is  the  elevation-angle  error. 


NUMERICAL  COMPARISON  OF  ONE-  AND  TWO-DIMENSIONAL  SOLUTIONS 
FOR  RECEIVERS  ABOVE  A  PLANE  EARTH 

The  one-dimensional  equations,  much  less  complicated  and  costly 
to  solve  numerically  than  the  two-dimensional  forms,  should  be  used 
whenever  possible.  Below,  we  quantify  the  altitude  regime  over  which 
the  simpler  equations  may  be  used  by  comparing  one-dimensional  solu¬ 
tions  with  those  of  the  more  accural?  two-dimensional  version.  To 
isolate  inaccuracies  due  solely  to  elevation,  we  make  the  comparisons 
for  a  plane  earth.  Section  IV  presents  a  complementary  analysis  that 
compares  one-  and  two-dimensional  solutions  for  terrain  irregularities 
of  various  sizes. 

Figure  1  shows  the  phase  of  W  calculated  from  the  one-  and  two- 

dimensional  equations  for  a  frequency  of  100  kHz  and  a  plane  earth 

_2 

having  a  conductivity  of  10  mhos/m.  The  curves  trace  phase  as  a 
function  of  distance  for  several  receiver  altitudes  between  1  m  and 
10  km. 


For  comparison,  we  also  show  the  phase  of  W  as  computed  from 
Norton's  ( 1 937 |  well-known  solution  of  the  differential  wave  equation — 
which  is  inherently  two-dimensional.  Although  Norton's  solution  is  a 
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tar  simpler  way  of  obtaining  W  for  a  plane  earth  than  is  either  of 
the  integral  equations,  it  is  valid  only  for  a  plane  earth,  whereas 
the  integral  equations  are  also  valid — subject  to  the  constraints 
given  in  Sec.  IV — for  propagation  over  irregular  terrain.  The  com¬ 
parison  in  Fig.  1  is  useful  for  defining  the  elevations  at  which  each 
integral  solution  agrees  with  ttie  classic  results  for  the  idealized 
model.  The  integral  equations  are  presumably  also  valid  at  the  same 
elevations  for  more  realistic  terrain,  where  Norton's  solution  cannot 
be  used. 

As  expected,  our  two-dimensional  solutions  agree  well  with 
Norton's  for  all  distances  and  elevations  shown.  The  small  differ¬ 
ences  are  due  to  rounding  in  our  numerical  solution  and,  at  very  short 

2  2 

distances,  to  Norton's  ignoring  terms  of  order  l/(kx^)  and  (z q/x^)  • 

Also  as  expected,  our  one-dimensional  results  are  accurate  provided 

the  constraints  on  z(.  expressed  by  Eq.  (20)  are  satisfied.  However, 

they  exhibit  large  errors  if  the  elevation  is  1  km  or  less  (end-point 

error)  or  if  7-q/xq  is  too  large  (elevation-angle  error).  The  error 

in  the  one-dimensional  results  for  z  =  1  m  is  so  great  as  to  be  off 

-1/2  U 

scale — a  consequence  of  the  z^  divergence  discussed  in  connection 
with  Eq.  (17). 

The  reliability  of  phase  predictions  made  using  the  one-dimen¬ 
sional  equation  is  illustrated  in  Fig.  2,  which  shows  contours  of 
constant  one-dimensional  errors  in  the  elevation-distance  plane.  We 
define  phase  errors  as  the  differences  (in  meters)  between  the  one- 
and  two-dimensional  results  shown  in  Fig.  1.  Figure  2  also  shows 
several  contours  of  constant  elevation  angle  z^/x^,  and  on  the  right 
axis,  the  elevation  z^  expressed  in  units  of  reduced  wavelength  X / 2tt . 
The  upper  branches  of  the  contours  correspond  to  elevation-angle 
errors;  the  lower  branches  correspond  to  end-point  errors. 

The  contours  in  Fig.  2  conform  closely  to  the  locus  of  points 
in  the  x^  -  z()  domain  corresponding  to  the  validity  conditions  of 
Eq .  (20).  For  example,  to  achieve  a  phase  error  smaller  than,  say, 

TO  m,  the  receiver  must  be  inside  the  50  m  contour.  From  Fig.  2,  we 
see  that  the  receiver  must  then  be  below  the  line  zq/xq  **  0.25 
(elevation-angle  error)  and  above  the  line  2ttXq/Zq  2  (end-point 
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error).  An  accuracy  of  25  m  requires  the  receiver  to  be  below  the 
line  zq/Xq  «  0.13  and  above  the  line  2ttZq/A  3.  The  usable  domain 
of  the  one-dimensional  equation  thus  shrinks  as  error  tolerance 
tightens. 

Fortunately,  the  end-point  and  elevation-angle  errors  caused  by 
the  one-dimensional  approximation  are  small  at  most  ranges  and  alti¬ 
tudes  relevant  to  airborne  receivers.  Johler  and  Berry  [1967],  for 
example,  present  one-dimensional  calculations  at  100  kHz  for  air¬ 
borne  receivers  at  distances  of  150  to  400  km  and  altitudes  of  2.5 
to  10  km;  Fig.  2  shows  that  range  of  parameters  to  be  almost  entirely 
within  the  10  m  error  contour.  However,  the  one-dimensional  equation 
can  yield  significant  phase  errors  at  100  kHz  for  receivers  positioned 
below  1  km  or — for  higher  altitudes — at  ranges  of  less  than  several 
tens  of  kilometers. 

In  practice,  much  of  the  error  for  receivers  that  are  above  the 
ground,  but  no  higher  than  1  or  2  km,  can  be  avoided — without  incurring 
the  expense  of  the  two-dimensional  solution — by  using  the  one-dimen¬ 
sional  formulation  for  grounded  receivers  even  though  the  receivers  are 
actually  slightly  elevated.  The  situation  is  illustrated  in  Fig.  3, 

which  again  pertains  to  a  frequency  of  100  kHz  and  a  plane  earth  of 
_2 

conductivity  10  mhos/m.  The  curves  shown  apply  to  all  ranges,  pro¬ 
vided  zq/xq  <  0.1  (to  make  elevation-angle  errors  unimportant).  For 
an  airborne  receiver  at  an  altitude  of  10  km,  for  example,  the  curves 
apply  to  all  ranges  beyond  100  km.  The  end-point  error  decreases 
with  elevation,  whereas  the  difference  between  the  phase  on  the  ground 
(computed  from  Eq .  (4) j  and  that  at  [computed  from  Eq .  (3)]  in¬ 
creases.  The'  two  curves  cross  at  1.3' km,  where  the  phase  error  is 
about  30  m. 

Thus,  at  least  for  the  conductivity  and  frequency  in  question, 

< ne-d imens iona 1  end-point  errors  can  be  kept  below  30  m  simply  by 
computing  the  phase  as  if  it  were  on  the  ground  [Eq.  (4)]  for  alti¬ 
tudes  under  1.3  km,  and  by  computing  the  phase  from  Eq .  (5)  for  greater 
heights.  The  accuracy  that  results  is  contained  within  the  shaded 
region  in  Fig.  3;  greater  accuracy  requires  numerical  solution  of 
the  two-dimensional  equation. 
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IV.  TWO-DIMENSIONAL  SOLUTIONS  FOR  ISOLATED  IRREGULARITIES 


We  showed  in  Sec.  Ill  and  Appendix  A  Chat,  even  if  the  terminals 
are  on  the  ground,  the  stationary-phase  derivation  of  the  one¬ 
dimensional  equation  is  valid  only  for  terrain  that  is  relatively 
uniform  across  the  width  of  the  first  Fresnel  zone,  and  that  more 
severe  irregularities  require  the  two-dimensional  equation.  For 
typical  pathlengths  of  a  few  hundred  kilometers,  the  Fresnel  zone 

1  /  7 

width  (Ax^)  is  on  the  order  of  tens  of  kilometers  at  a  frequency 
of  100  kHz.  Since  irregularities  in  actual  terrain  are  often  smaller 
than  tens  of  kilometers,  the  regime  of  the  one-dimensional  equation's 
applicability  is  quite  restricted. 

This  section  compares  numerical  solutions  of  the  one-  and  two- 
dimensional  equations  for  propagation  over  surfaces  having  irregular¬ 
ities  smaller  than  a  Fresnel  zone.  First,  we  modify  the  two-dimen¬ 
sional  equation  to  a  form  more  efficient  for  computing  the  effects  of 
isolated  irregulari ties  than  the  standard  form,  Eq .  (2).  Then  we 
give  an  approximate  analytic  solution  to  this  equation  which — although 
valid  only  under  limited  conditions  —  lends  insight  to  the  dependence 
of  W  on  the  location  and  transverse  dimension  of  the  irregularity. 
Finally,  we  present  solutions  to  the  full  two-dimensional  equation. 

INTEGRAL  EQUATIONS  FOR  ISOLATED  IRREGULARITIES  ON 
A  HOMOGENEOUS  PLANE 

The  standard  form  of  the  two-dimensional  integral  equation — given 
by  Eq.  (2) — requires  integration  over  the  entire  surface.  Such  an 
integration  is  necessary  for  continuously  varying  or  undulating  ter¬ 
rain,  hut  not  for  terrain  that  is  uniform  except  for  a  few  isolated 
irregularities.  In  the  case  of  mostly  regular  terrain,  we  can  take 
auvantage  of  the  fact  that  the  attenuation  function  is  already  known 
on  a  plane  or  sphere  having  the  properties  of  the  homogeneous  part 
o!  the  surface.  Specifically,  we  write 


W  =  w}wH 


(21) 
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where  WH  is  the  known  solution  for  a  uniform  surface,  and  W  is  the 

new  unknown  function.  For  propagation  on  a  spherical  earth,  W  woulc 

H 

be  the  well-known  residue-series  solution  [Bremmer,  1958].  Here,  we 

consider  irregularities  on  an  otherwise  uniform  plane  for  which  W  is 

H 

the  Norton-Sommerf eld  function  {Norton,  1937]  given  by 


W^(r)  =  1  +  i  /irn  e  ^  erfc(-i  i/q)  , 


(22) 


where 


Tl  =  ikr /2n^  (23) 

and  erfc  is  the  complement  of  the  error  function  [Abromowitz  and 
Stegun,  1972]. 

Consider  first  a  plane  where  the  refractive  index  is  n  every¬ 
where  except  for  a  certain  region  where  it  is  n^ — a  lake  on  flat  land 
or  a  flat  island  in  the  ocean,  for  example.  In  such  a  case.  Appendix  B 
shows  that  the  integral  equation  governing  is 


Vro>  ' 


1  + 


ik 

2tt 


dx  dy  A(x,  y) 


ik(r1+r2-r()) 


rlr2 


X 


~“H(ri)HH<r2) 

WH(r0) 


Wj(x,  y)  , 


(2A) 


where 


A(x,  y)  =  - - - -  (25) 

ng(x,  y)  g 

is  the  impedance  contrast  between  the  lake  or  island  and  the  surround¬ 
ing  region. 

At  first,  Eq.  (24)  looks  more  formidable  than  Eq.  (2)  because 
W  appears  in  the  integrand.  It  is,  however,  much  less  costly  to 
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solve  numerically  because  the  integration  need  be  carried  out  only 
over  regions  where  A  ^  0,  rather  than  over  the  entire  x  -  y  plane. 

The  savings  are  well  worth  the  added  complexity  of  including  the 
Norton-Sommerfeld  functions  in  the  integrand.  Equation  (24)  is  re¬ 
lated  to  a  two-dimensional  integral  equation  derived — but  not  solved — 
by  Feinberg  [1945]. 

As  a  second  example,  consider  a  uniform  plane  of  refractive  in¬ 
dex  n^  on  which  sits  an  isolated  hill.  Unlike  the  first  example,  we 
have  no  impedance  contrast  here,  but  have  instead  a  region  where 
3r2/3n  t  0.  For  such  a  case,  we  show  in  Appendix  B  that  satisfies 


VV  =  1  // dx  dy 


ik(r1+r2-r0) 


WH(rl)WH(r2) 

VV 


Wx(x,  y)  , 


(26) 


which  is  similar  to  Eq.  (24)  except  that  terms  involving  dr^/on  appear 
in  place  of  terms  involving  A.  Equation  (26)  is  numerically  more  con¬ 
venient  than  Eq.  (2)  because  integration  is  required  only  where 
<3r2/3n  is  significant.  However,  the  integration  region  can  extend 
well  beyond  the  hill  itself;  Eq .  (26)  is  thus  more  costly  to  solve 
than  Eq .  (24),  which  requires  integration  only  over  the  region  actually 
occupied  by  the  irregularity.  That  is,  the  "irregularity"  caused  by 
an  isolated  hill  is  characterized  bv  3r2/3n,  which  affects  a  region 
much  larger  than  the  hill  itself. 


BORN-APPROXIMATION  SOLUTION  FOR  WEAK,  ISOLATED  IRREGULARITY 

We  can  solve  Eq .  (24)  approximately  for  a  weak  impedance  irregu¬ 
larity  embedded  in  a  perfectly  conducting  plane,  where  is  unity 
and  W ^  =  W.  We  set  «  1  on  the  right  side  to  obtain 


W,  «  1  + 


ik 

2n 


// 


dx  dy  A(x,  y) 


ik(r1+r2-xQ) 


rlr2 


(27) 
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for  terminals  on  the  x  axis.  Equation  (27)  is  valid  if  the  second 
term  on  the  right  is  small,  giving  a  posteriori  justification  of  our 
assumption  that  a*  1.  If  the  irregularity  were  a  small  hill  not 
too  near  the  terminals,  3r^/3n  would  replace  A  in  Eq.  (27).  This 
perturbation  solution  is  equivalent  to  the  well-known  Born  approx¬ 
imation. 

Equation  (27)  requires  a  double  integration  over  the  impedance 
contrast  A(x,  y) .  We  assume  the  simple  form 

MX,  ,)  -  M  «P  [-  -  ^iiil  .  (28) 

U  L  (Axr  (Ay)  J 

which  allows  us  to  illustrate  the  dependence  of  on  the  lateral 
gradients  and  position  of  the  perturbation  by  adjusting  (Ax),  (Ay), 
x,  or  y,  respectively.  We  can  perform  the  integrations  using  the 
saddle-point  approximation,  provided  kx^  »  1,  Ax  «  x^,  and 
2Ax  <  x  <  2(Xq  -  Ax).  The  procedure  is  described  by  Field  and 
Joiner  (19791,  and  the  result  is: 


W1  w  1  +  iAQ(Ax) 


ikx„ 


2x(Xq  -  x) 


1/2 
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(29) 


where 


and 


A(x0,  x.  Ay) 
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1  +  i  (  a\2 

4x(xfl  -  x)] 

it  \Ay  / 

- 1 

CM  O 

X 

_ 1 

d 


(30) 


(31) 


In  Eqs.  (29)  through  (31),  d  is  the  maximum  half-width  of  the  first 
Fresnel  zone.  A  accounts  for  the  width  and  longitudinal  position  of 


X  - 
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the  perturbation,  and  the  exponent  in  Eq.  (29)  for  its  transverse 

—  —  2 

location.  The  factor  4x(x^  -  x)/Xq  in  Eq.  (30)  accounts  for  the 
narrowing  of  the  Fresnel  zone  away  from  midpath,  and  has  a  maximum 
value  of  unity  when  x  =  x^/ 2.  Although  valid  only  under  very  re¬ 
stricted  conditions,  Eqs.  (29)  through  (31)  reveal  several  important 
characteristics  of  the  attenuation  function. 

Weak  One-Dimensional  Irregularity 

By  taking  the  limits  Ay  »  d  and  y,  we  find  that  Eq.  (29)  assumes 
the  following  form  for  a  wide  (i.e.,  one-dimensional)  on-path  pertur¬ 
bation  : 


1  +  A 


(Ax)  (ikx0)1,/2 


Av-*3" 


°  [2x(x0  '  x)] 


-^1/2  • 


(32) 


The  denominator  in  Eq.  (32)  is  largest  at  midpath,  x  =  x„/2,  and 

*  —  u 

small  near  the  endpoints,  x  «  0  or  xQ.  The  effects  of  the  perturba¬ 
tion  are  therefore  greatest  if  it  is  near  either  terminal,  and 
smallest  near  midpath.  This  enhanced  importance  of  regions  near  the 
terminals  is  well  known  [Feinberg,  1959]. 

Narrow  On-Path  Irregularity 

For  perturbations  much  narrower  than  a  Fresnel  zone  Ay  <<  d,  and 
Eq .  (30)  becomes 

A  _  „.iii(.AX/d)2  _  (33) 

Ay<<d  4x(Xq  -  x)/Xq 

The  second  term  on  the  right  of  Eq.  (32)  must  therefore  be  reduced 
by  a  factor  of  magnitude 

* 

Recall  that  the  derivation  of  Eq.  (29)  prohibits  moving  the 
perturbation  to  within  ~2Ax  of  either  terminal. 


-25- 


2yx(xQ  -  x)/Xq 


(34) 


to  account  for  the  limited  width  of  the  perturbation.  Equation  (34) 
gives  the  intuitively  reasonable  result  that  the  two-dimensional 
correction  factor  is  proportional  to  the  fraction  of  the  Fresnel  zone 
effectively  occupied  by  the  impedance  contrast. 

For  narrow  perturbations  the  reduction  indicated  by  Eq.  (34) — 
but  omitted  by  one-dimensional  equations — can  be  considerable.  For 
example,  for  Ay  =  3  1cm  (a  fairly  sizable  terrain  feature),  X  =  3  km, 

Xq  =  500  km,  and  x  =  xq/2,  this  two-dimensional  correction  factor  is 
about  0.27.  For  these  parameters,  the  one-dimensional  formulation 
therefore  overstates  the  impact  of  the  perturbation  by  nearly  a  fac¬ 
tor  of  four. 

Although  strictly  valid  only  for  the  weak  Gaussian  impedance 
contrast  assumed  above,  the  above  discussion  shows  that  the  one¬ 
dimensional  integral  equation  is  in  serious  error  unless  the  pertur¬ 
bation  is  nearly  as  wide  as  a  Fresnel  zone.  Details  will  depend 
on  the  structure  of  a  given  irregularity,  but  the  conclusion  is  gen¬ 
eral;  viz.  two-dimensional  solutions  are  needed  unless  terrain  fea¬ 
tures  are  nearly  uniform  over  transverse  distances  on  the  order  of 
/AXq.  The  notion  that  the  stationary-phase  reduction  of  the  two- 
dimensional  equations  requires  slight  variation  across  a  Fresnel  zone 
is  certainly  not  new  [Bremmer,  1958;  Feinberg,  1945].  Nonetheless, 
the  inaccuracy  of  the  one-dimensional  equation  at  low  frequencies  for 
terrain  features  narrower  than,  say,  ten  kilometers  is  often  overlooked 
in  applications. 

Next,  we  examine  the  dependence  of  on  pathlength  by  rewriting 
Eq.  (34)  in  the  form 
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which  shows  that,  for  fixed  x,  the  effect  of  the  perturbation  diminishes 
with  increasing  pathlength  according  to  [ i  —  x/x^i  ^2.  This  "re¬ 
covery"  of  Wj  occurs  because  the  width  of  the  Fresnel  2one  increases 
as  the  path  is  made  longer;  and  the  fraction  occupied  by  the  perturba¬ 
tion  therefore  decreases.  This  two-dimensional  phenomena  causes  W 
to  recover  beyond  an  obstacle  more  quickly  than  predicted  by  the  one- 
dimensionaJ  approximation. 

Isolated  Off-Path  Irregularity 

We  showed  above  that  the  one-dimensional  result,  Eq.  (32),  must 
be  reduced  by  the  factor  in  Eq .  (34)  to  account  for  an  on-path  per¬ 
turbation  being  narrower  than  a  Fresnel  zone.  If  a  narrow  terrain 
feature  is  centered  off-path,  the  result  in  Eq.  (32)  must  be  reduced 
by  still  another  factor.  By  combining  Eq.  (30)  with  the  exponential 
term  in  Eq.  (29),  we  find  this  "off-path"  factor  F  to  be 


F  ss  exp 


-2 

X- 


(Ay)2  + 


(36) 


where,  for  simplicity,  we  have  assumed  x  =  Xq/2.  Equation  (36)  shows 
that  the  "off-path"  factor  is  governed  by  the  Fresnel  zone  width  for 
narrow  perturbations  (Ay  «  d),  and  depends  only  weakly  on  off-path 
distance  y,  for  wide  ones.  F  is  related  to  a  more  complicated  off- 
path  factor  used  by  King  and  Wait  [1976]. 


NUMERICAL  SOLUTIONS  OF  THE  TWO-DIMENSIONAL  EQUATION 

The  algorithm  used  to  solve  Eqs.  (24)  and  (26)  is  given  in 
Appendix  C.  The  first  two  examples  given  below  compare  the  numerical 
solutions  against  analytic  and  experimental  results,  thereby  validat¬ 
ing  our  approach.  The  third  example  compares  solutions  of  the  one- 
and  two-dimensional  equations  for  terrain  features  of  various  widths. 
Because  computational  expense  is  a  consideration,  all  examples  per¬ 
tain  to  impedance — rather  than  topographic — irregularities. 


-27- 


Weak  Gaussian  Impedance  Contrast 

Equations  (29)  and  (32)  give  two-  and  one-dimensional  analytic  so¬ 
lutions,  respectively,  for  a  weak  Gaussian  impedance  contrast  ( Kq .  (26)) 
in  a  conducting  plane.  Their  main  value  is  to  lend  insight  into  the 
behavior  of  W,  since  their  range  of  validity  is  too  restricted  for 
application  to  actual  terrain.  Nonetheless,  if  we  confine  ourselves 
to  weak  perturbations  where  |Wj  -  l|  <<  1,  we  can  use  Eqs.  (29)  and 
(32)  to  validate  numerical  solutions  of  Eqs.  (24)  and  (4),  respectively. 

The  following  table  gives  the  phase  of  as  computed  by  the 

analytic  approximations  and  numerical  solutions.  The  assumed  maximum 

_2 

impedance  contrast  A^  corresponds  to  a  conductivity  0  of  10  mhos/m, 
and  the  frequency  is  100  kHz.  The  other  parameters  (see  table)  were 
chosen  to  satisfy  the  validity  restrictions  on  Eqs.  (29)  and  (32). 

If  we  had  selected  parameters  to  give  a  more  substantial  — e.g., 

greater  Ax  or  lower  conductivity — we  would  have  violated  these  re¬ 
strictions.  Moreover,  since  */<Xx^/2  8.7  for  the  case  shown,  our  use 

of  Ay  =  1C  km  corresponds  to  a  relatively  wide  perturbation  that  does 
not  cause  major  two-dimensional  effects.  Use  of  a  much  narrower  per¬ 
turbation  would  reduce  -  1  to  where  it  could  not  be  resolved  within 
our  numerical  accuracy. 


COMPARISON  BETWEEN  ANALYTICAL  AND  NUMERICAL  SOLUTIONS 
FOR  A  WEAK  GAUSSIAN  IMPEDANCE  CONTRAST 


Phase  of 

(rad) 
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For  the  above  reasons,  the  phase  anomalies  shown  in  the  table 
are  very  small  and  the  one-  and  two-dimensional  results  do  not  differ 
substantially.  Nonetheless,  these  results  confirm  two  important  facts. 
First,  our  numerical  algorithms  are  accurate.  Second,  the  one-dimen¬ 
sional  equation  overstates  the  effect  of  an  on-path  irregularity  that 
is  narrower  than  a  Fresnel  zone. 

Conducting  Region  on  Dielectric  Plane 

King  and  Tsukamoto  [1966]  made  laboratory  measurements  of  surface 
waves  propagating  over  a  plane  dielectric  surface  containing  isolated 
conducting  regions.  De  Jong  [1975]  obtained  good  agreement  with  these 
data  by  numerically  solving  a  two-dimensional  integral  equation  similar 
to  our  Eq.  (24),  albeit  of  somewhat  different  form.  To  our  knowledge, 
his  results  are  the  only  full  two-dimensional  solutions  available  in 
the  literature.  Before  proceeding  to  models  representing  groundwave 
propagation,  we  compare  our  calculations  with  de  Jong’s. 

Figures  4  and  5  show  the  propagation  model  and  our  numerical  solu¬ 
tions  for  the  amplitude  and  phase  of  W^.  A  rectangular  "island"  of 

infinite  conductivity  is  embedded  in  a  dielectric  plane  having  an 
-1  -3 

impedance  (n  )  of  0.11  -  9.35  *  10  i.  The  width  and  length  of 
£ 

the  island  are  3  and  6.75  km,  respectively,  and  the  frequency  is 
100  Hz.  The  width  of  the  Fresnel  zone  is  /3  x  30  «  9.5  km  for  a  path- 
length  of  30  km.  The  island  is  therefore  much  narrower  than  a  Fresnel 
zone,  and  we  expect  large  differences  between  the  one-  and  two-dimen¬ 
sional  solutions. 

Our  two-dimensional  results  are  in  excellent  agreement  with 
de  Jong's,  if  allowance  is  made  for  hi’’  assumed  transmitter  and  re¬ 
ceiver  being  a  small  fraction  of  a  wavelength  above  the  surface. 
Moreover,  they  agree  well  with  the  data  (not  repeated  here)  every¬ 
where  except  in  a  narrow  region  beginning  at  the  rear  of  the  island 
and  extending  for  one  or  two  wavelengths.  Even  our  calculated  standing- 
wave  structure  agrees  closely  with  King  and  Tsukamoto 's  data  and 
de  Jong's  calculations. 

We  concur  with  de  Jong's  explanation  that  the  region  of  poor 
agreement  with  experiment  is  caused  by  the  well-known  failure  of 


Distance  (km) 


Fig.  ^-Calculated  magnitude  of  W]  for  a  conducting  "island"  (shaded 
region)  on  a  dielectric  plane:  frequency  =  100  kHz, 
nonelevated  terminals 
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impedance  boundary  conditions  near  abrupt  conductivity  changes.  The 
discontinuous  behavior  of  the  phase  at  the  rear  of  the  island  indi¬ 
cates  incorrect  boundary  conditions.  It  is  encouraging  that  this  poor 
accuracy  persists  for  only  one  or  two  wavelengths  beyond  the  boundary, 
after  which  our  results  (and  de  dong's)  recover  to  the  correct  values. 

Figures  4  and  5  quantify  several  defects  of  the  one-dimensional 
solution  which  are  expected  on  intuitive  grounds.  ’.irst,  the  one¬ 
dimensional  solution  overstates  the  amplitude  anomaly  caused  by  the 
island — a  consequence  of  its  treating  the  island  as  infinitely  wide. 
Second,  it  recovers  too  slowly  beyond  the  island — a  consequence  of 
not  accounting  for  Fresnel  zone  widening  as  the  pathlength  increases. 
Third,  it  omits  the  standing-wave  patterns  evident  from  the  two- 
dimensional  solutions — a  consequence  of  the  stationary-phase  approx¬ 
imation  neglecting  terrain  features  beyond  the  receiver. 

Recall  that  Figs.  4  and  5  pertain  to  W^,  which  accounts  solely 
for  the  effect  of  the  island.  These  results  must  be  multiplied  by 
W  [see  Eq.  (21)],  given  in  Fig.  6,  to  obtain  the  full  attenuation 

n 

function  W. 

Saltwater  "Lake"  on  Poor  Ground 

Having  verified  our  approach  for  two  idealized  cases,  we  next 
consider  parameters  more  representative  of  long-wave  propagation  on 
the  earth.  Even  here  we  use  a  simple  model  to  avoid  unnecessary  com¬ 
puter  expense.  Our  method  can  readily  treat  more  complicated  situa¬ 
tions,  albeit  at  increased  cost. 

Our  iliustrative  model  assumes  a  conductivity  <i(x,  y)  given  by 
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Fig.  6--Attenuation  function  for  plane  earth:  frequency  =  100  kHz 
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This  conductivity  might  correspond  to  a  saltwater  lake  or  bay  of  vari¬ 
able  depth.  We  assume  this  "lake"  to  be  embedded  in  ground  of  con- 
-3 

ductivitv  10  mhos/m  and  dielectric  constant  10,  and  to  extend  from 
x  =  40  km  to  x  =  60  km.  Its  half-width  is  specified  bv  Ay,  which  we 
will  vary.  The  model  given  by  Eq .  (37)  has  the  desirable  feature  of 
continuous  lateral  derivatives.  However,  our  assumed  .arge  conduc¬ 
tivity  contrast  causes  significant  impedance  changes  within  a  fraction 
of  a  kilometer  of  the  boundaries. 

Figure  7  shows  two-dimensional  solutions  |  F.q .  (24)  |  for  W^  versus 
distance  for  half-widths  Ay  of  3,  5,  and  7  km.  We  also  show  for  refer¬ 
ence  the  one-dimensionai  solution  | Eq .  (4)].  All  show  little  or  no 
effect  in  front  of  the  lake,  a  pronounced  enhancement  in  it,  and  various 
degrees  of  recovery  beyond  it.  The  levels  of  enhancement  and  recovery 
shown  should  >u:t.  be  generalized  to  other  parameters.  Wait  [1964]  gives 
comprehensive  one-dimensional  results  for  land-sea-land  propagation, 
and  shows  that  the  enhancement  in  the  sea  portion,  and  subsequent  re¬ 
covery,  depend  strongly  on  virtually  all  parameters  of  the  problem. 

Even  simple-looking  one-dimensional  results  shown  in  Fig.  7  could 
change  drastically  if  we  altered  the  lake's  location  or  length,  or  the 
ground's  conductivity.  The  results  in  Fig.  7  must  be  multiplied  by 

W  ,  given  in  Fig.  6,  to  recover  the  full  attenuation  function  W. 

H 

The  accuracy  of  the  solutions  is  poor  near  the  rear  boundary, 
because  the  impedance  boundary  conditions  fail  near  strong  lateral 
gradients.  However,  the  preceding  example  showed  that  the  accuracy 
recovers  at  some  distance  beyond  the  irregularity.  Thus,  we  should 
restrict  attention  to  regions  at  least  a  few  wavelengths  from  the 
boundaries  of  our  model  lake. 

The  maximum  half-width  of  the  first  Fresnel  zone  is  about  12 
and  17  km  for  pathlengths  of  200  and  400  km,  respectively.  All  models 
shown  in  Fig.  7  are  therefore  fairly  narrow,  although  the  one  having 
Ay  =  7  km  occupies  a  sizable  fraction  of  the  first  zone.  As  expected, 
the  two-dimensional  solutions  more  closely  approach  the  one-dimensional 
solution  as  Ay  increases.  Also,  as  predicted  by  Eq .  (35),  the  one¬ 
dimensional  solution  recovers  too  slowly  beyond  the  lake.  Finally, 
and  most  important,  the  one-dimensional  solution  erroneously  gives 
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-Magnitude  of  W]  for  seawater  regions  (shaded  area)  of  various  widths 
on  a  poorly  conducting  plane  earth:  frequency  =  100  kHz 
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the  same  result  for  all  values  of  Ay.  Even  for  Ay  =  7  km — wider  than 
many  natural  impedance  or  topographic  Irregularities — the  one-dimen¬ 
sional  solution  strongly  overstates  the  propagation  anomaly  for  long 
pathlengths . 

Figure  7  shows  that,  although  the  one-  and  two-dimensional  solu¬ 
tions  approach  one  another  for  wide  disturbances  and  long  pathlengths, 
they  differ  for  several  tens  of  kilometers  beyond  the  lake.  Thus,  we 
have  the  apparent  paradox  that — near  a  terrain  feature — the  solution 

of  the  two-dimensional  equation  dour,  not  approach  that  of  the  one- 

* 

dimensional  equation  as  Ay  ■>  °°.  We  believe  that  this  disagreement 
occurs  because  the  two-d imensional  equation  is  more  accurate,  even  for 
a  one-dimensional  terrain  feature — e.g.,  one  that  has  no  y-dependence. 
One-dimensional  features  cause  reflections,  resonances,  and  standing 
waves — all  of  which  are  lost  in  the  stationary-phase  approximation,  but 
retained  in  the  two-dimensional  formulation.  Such  phenomena  diminish 
in  importance  as  the  distance  from  boundaries  becomes  large,  and  the 
one-dimensional  result  should  be  accurate  for  large  values  of 
Xq  -  (x  +  Ax).  This  behavior  is  evident  in  Fig.  7.  Of  course,  the 
accuracy  of  even  the  two-dimensional  equation  is  adversely  affected  by 
the  failure  of  impedance  boundary  conditions  near  abrupt  conductivity 
changes . 

Figures  8  and  9  show  as  a  function  of  transverse  receiver 
position  y  for  several  values  of  x  and  lakes  having,  respectively, 
half-widths  of  3  and  7  km.  The  function  exhibits  a  classic  dif¬ 
fraction  pattern  of  maxima  and  minima  in  the  transverse  direction. 

That  pattern  is  to  be  expected  because  the  obstacle  is  smaller  than 
a  Fresnel  zone.  Although  the  details  of  the  signal  structure  shown 
in  Figs.  7  through  9  pertain  to  the  model  and  geometric  factors  assumed, 
a  qualitatively  similar  structure  will  occur  for  any  terrain  feature 
having  lateral  dimensions  smaller  than  a  Fresnel  zone.  Diffraction 
patterns  such  as  those  shown  in  Figs.  8  and  9  are,  of  course,  omitted 
by  the  widely  used  one-dimensional  equation. 

This  conclusion  does  not  contradict  Eq.  (32)  and  the  results  in 
the  table  (see  p.  27),  which  are  valid  only  for  perturbations  too 
weak  to  cause  significant  reflection. 
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V.  CONCLUSIONS 


The  two-dimensional  integral  equation  for  groundwave  propaga¬ 
tion  has  different  forms  for  elevated  and  ground-based  receivers. 

The  two  forms  become  equivalent  in  the  limit  of  zero  elevation.  The 
stationary-phase  approximation  can  sometimes  be  used  to  derive  much 
simpler  one-dimensional  equations,  which  are  widely  used  for  most 
applications.  This  stationary-phase  reduction  in  dimensionality  causes 
so-called  elevation-angle  and  end-point  errors,  both  of  which  are 
associated  with  elevation.  We  have  quantified  these  errors  by  com¬ 
paring  our  one-  and  two-dimensional  numerical  solutions  with  Norton's 
attenuation  function  for  a  uniform  flat  earth.  The  two-dimensional 
solutions  agree  closely,  whereas  the  one-dimensional  solutions  fail 
at  certain  elevations. 

F.levat ion-angle  errors  are  large  unless  the  elevation  is  much 
smaller  than  the  pathlength.  End-point  errors  occur  because  the 
stationary-phase  approximation  improperly  accounts  for  the  region 
near  the  receiver,  and  are  large  for  elevated  receivers  at  heights 
below  about  one-sixth  of  a  wavelength.  The  one-dimensional  approxi¬ 
mations  for  grounded  and  elevated  receivers  are  therefore  discontinuous 
near  the  ground.  Fortunately,  neither  type  of  error  is  too  serious  at 
a  frequency  of  100  kHz  for  most  ranges  and  altitudes  pertaining  to 
airborne  receivers.  However,  end-point  errors  are  significant  at 
altitudes  below  1  to  2  km,  and  elevation-angle  errors  can  be  important 
for  pathlengths  of  tens  of  kilometers  or  less.  In  practice,  for  re¬ 
ceivers  that  are  above  the  ground  but  no  higher  than  1  or  2  km,  much 
of  the  end-point  error  can  be  avoided — without  the  expense  of  the  two- 
dimensional  solution — simply  by  using  the  one-dimensional  solution 
for  grounded  receivers. 

Our  numerical  solutions  of  the  two-dimensional  equation  for  iso¬ 
lated  terrain  features  agree  well  with  (1)  approximate  solutions  from 
first-order  perturbation  theory,  (2)  numerical  solutions  of  de  Jong 
[19751  for  an  idealized  model,  and  (3)  laboratory  measurements  by 
King  and  Tsukamoto  [1966]. 
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It  is  well  known  that  the  one-dimensional  equation  is  invalid 
unless  the  terrain  is  nearly  uniform  across  a  Fresnel  zone.  Compari¬ 
sons  between  our  one-  and  two-dimensional  solutions  for  on-path 
terrain  quantify  the  errors.  For  a  frequency  of  100  kHz  and  features 
narrower  than  about  10  km,  we  find  that  the  one-dimensional  equation 
erroneously  predicts  propagation  anomalies  that  (1)  are  independent  of 
width,  and  therefore  too  large,  (2)  diminish  too  slowly  at  long  dis¬ 
tances,  and  (3)  do  not  exhibit  a  diffraction  pattern  in  the  transverse 
direction. 

The  one-  and  two-dimensional  solutions  approach  one  another  at 
large  distances  beyond  very  wide  terrain  features.  However,  they 
disagree  near  boundaries  because  the  stationary-phase  approximation 
neglects  reflections  and  interference  effects.  Therefore,  even  for 
terrain  that  exhibits  no  transverse  gradients,  only  the  two-dimen¬ 
sional  equation  accounts  for  detailed  signal  structure. 

Considerable  error  can  be  incurred  at  low  frequencies  by  apply¬ 
ing  the  one-dimensional  equation  to  moderately  sized  terrain  features. 
For  example,  for  a  pathlength  of  500  km,  that  equation  overstates  by 
a  factor  of  4  the  effect  of  an  obstacle  6  km  in  diameter.  It  cannot 
give  accurate  results  unless  the  diameter  approaches  a  Fresnel  zone 
width,  which  for  this  example  is  several  tens  of  kilometers.  None¬ 
theless,  the  one-dimensional  equation  is  often  applied  to  terrain  that 
does  not  satisfy  the  validity  requirements. 

For  the  above  reasons,  we  believe  the  one-dimensional  equation 
to  be  incapable  of  treating  many  terrains,  even  if  the  input  data 
were  perfect.  Unfortunately,  routine  application  of  the  more  accurate 
two-dimensional  equation  to  irregular  terrain  is  probably  impractical. 
A  more  frui: ^ul  approach  is  to  devise  means  of  averaging  terrain  over 
the  Fresnel  zone  to  obtain  equivalent  one-dimensional  models.  The 
statistical  approach  of  Feinberg  [1944]  can  be  used  to  derive  a  one¬ 
dimensional  ;uation  for  the  average  field  if  the  terrain  has  many 
small,  randomly  located  irregularities.  His  approach  is  invalid  for 
terrain  exhibiting,  say,  a  few  hills  or  lakes,  however. 
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Appendix  A 

STATIONARY-PHASE  APPROXIMATION  FOR  ELEVATED  RECEIVERS 


For  terminals  on  a  plane  earth,  the  phase  of  the  integrand  in 
Eq .  (2)  is  dominated  by  the  term  exp  [ik(r^  +  r0  -  r^) ] ,  whose  phase 
is  constant  on  a  family  of  ellipses  defined  by  the  expression  r^  + 
constant,  but  varies  rapidly  along  the  hyperbolas  normal  to  the 
ellipses  (see  Fig.  A.l).  That  behavior  led  Hufford  to  use  elliptical 
coordinates  in  his  solution;  the  hyperbolas  shown  in  Fig.  A.l  define 
his  transverse  coordinate.  The  main  contributions  to  the  integral 
occur  within  the  ellipse  defined  by  r^  +  x^  =  A/2 — that  is,  within 
the  first  Fresnel  zone  (illustrated  in  Fig.  A. 2). 

Ellipses  and  hyperbolas  are  natural  coordinates  for  a  receiver 
on  a  plane  earth,  but  not  for  an  elevated  receiver.  For  the  elevated 
case,  we  must  begin  with  Eq .  (3),  which  contains  the  term  exp  [ik(r^  + 
R2  ~  V1  3n<*  ^°r  t*le  relation  rx  +  R2  =  constant  does  not  trace 

ellipses  in  the  plane  x  =  0.  Another  objection  to  elliptical  coordi¬ 
nates  is  that  for  either  ground-based  or  elevated  receivers,  they  are 
awkward  for  computing  higher  order  corrections  to  the  stationary-phase 
integration. 

Field  and  Allen  [1978],  in  rederiving  Eq.  (2)  for  nonelevated 
receivers,  used  cartesian  coordinates  and  applied  the  stationary-phase 
approximation  along  the  y  coordinate.  Despite  those  differences  from 
Hufford 's  approach,  the  two  results  agree.  The  reason  is  that  for 
either  method  to  be  valid  in  the  case  of  ground-based  receivers, 
kr^  »  1  must  hold,  which  in  turn  implies  r^  »  /Ar^.  The  last  con¬ 
dition  stretches  the  Fresnel  zone  into  a  highly  elongate.'  ellipse. 

The  hyperbolic  trajectories,  except  near  the  ends,  then  conform  nearly 
to  the  paths  defined  by  x  =  constant.  Hence,  integrating  along  the 
hyperbolas  is  nearly  equivalent  to  integrating  along  the  y  coordinate 
over  most  of  the  path. 

STATIONARY-PHASE  RESULT  FOR  ELEVATED  RECEIVERS 

To  assess  the  end-point  and  elevation-angle  errors  in  the  one¬ 
dimensional  equation,  we  extend  the  stationary-phase  integration  in 


1  N 

!  Hyperbolas  along  which  r  1  +  r2 

change  most  rapidly. 

Fig.  A. 1--E1 liptical  coordinate  system 
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cartesian  coordinates  so  that  it  applies  to  elevated  receivers,  then 
evaluate  correction  terms  to  the  usual  result.  By  performing  the 
derivation  for  a  plane  earth,  we  concentrate  solely  on  effects  related 
to  receiver  elevation.  We  seek  to  determine  when  the  two-dimensional 
equation  for  elevated  receivers. 


2W(P)  = 


+ 


ik(ri+R2-R0) 

kR2/R2  6 


(9) 


can  be  approximated  by  the  usual  one-dimensional  formula 


2W(x 


0’ 


V 


dx 


xR^(x  +  R2> 


1/2 


x  W(x,  0) 


+ 


ik(x+R2-RQ) 

e 


(A.l) 


Equation  (9)  is  simply  repeated  from  the  text;  F.q.  (A.l)  results 
from  inserting  Eqs.  (6)  and  (8)  into  Eq.  (5)  and  using  a  cartesian 
coordinate  system  with  the  ground  plane  at  z  =  0,  the  transmitter  at 
the  origin,  and  the  receiver  at  x^,  For  Eq.  (9),  the  integration 

point  Q  has  coordinates  x,  y,  0,  and  the  following  re 1: t ions  hold: 


and 


(A. 2) 


2 

y  ; 


(A. 3) 


x)2  + 


(A. 4) 


2  2 
r2  +  z0 


(A. 5) 
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For  the  one-dimensional  equations,  the  same  relations  apply,  with  y 
set  equal  to  zero. 

Equation  (9)  can  be  written 


°o  oo  ^  ^  fc. 

f  dx  f  dy  r (x,  y)  exp  |  ikR0  — ~  ~  ^ 

•'-oo  •'-oo  L  ® 


2W(Xq,  zq)  =  1  + 


(A.  6) 


where 


r(x 


•  y)  •  2i(vO"<X'  y'  0>  ["1  '  i1  +  k0Rj] 


(A. 7) 


The  inner  integral  in  Eq.  (A. 6)  is  given  by 


Kx) 


OO 

L 


dx  T(x,  v)  e 


ikRph(x ,y) 


(A. 8) 


where 


h  = 


rl  + 


1  . 


(A. 9) 


For  long  transmission  paths,  kr^  >>  1,  and  the  integral  [F,q.  (A. 8)} 
is  of  the  classic  form  amenable  to  approximate  evaluation  by  the 
stat ionary-pliase  method  (Erdelyi,  1956],  provided  T(x,  y)  varies  more 
slowly  than  exp  (ikKgh).  If  it  does,  the  stationary-phase  formula 
can  be  used  to  write 


I  (x)  e 


TTi/4 


ikR0h  (_22lY/2 

\kR()h"/ 


y=0 


(A. 10) 


By  inserting  Eq.  (A. 10)  and  the  relations 
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i 


and 


R0h 


=  X  +  R„ 


y=0 


(A. 11) 


(A. 12) 


into  Eq.  (A. 6),  we  recover  the  standard  form  of  the  one-dimensional 
equation  [Eq.  (A.l)].  The  x  integration  is  terminated  at  0  and  xQ, 
rather  than  at  ±°°,  because — as  discussed  below — the  exponent  varies 
slowly  for  0  £  x  £  x^  but  oscillates  rapidly  for  x  <  0  or  x  >  x^. 


END-POINT  ERROR  IN  ONE-DIMENSIONAL  APPROXIMATION 

The  above  treatment  is  valid  provided  T(x,  y)  varies  only  slightly 
over  distances  the  order  of 


<5y 


(kR0h") 


-1/2 


y=0 


(A. 13) 


That  condition  is  satisfied  for  most  values  of  x.  However,  nt  the 
receiver,  where  x  x^  and  R£  «  zQ,  Eq.  (A. 13)  gives 


<5y 


(A. 14) 


as  the  distance  over  which  T  must  be  relatively  constant;  whereas  the 
term 


i  !0 
kR2  1*2 


k(r2  +  zo) 


(A. 15) 


which  appears  in  T  [see  Eq.  (A. 7)],  has  a  peak  at  x  =  x^  and  a  half¬ 
width  of  Zg.  Therefore,  unless 
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z 


0 


or,  equivalently. 


Zg  >  A/2tt  ,  (A. 16) 

the  stationary-phase  approximation  fails  near  the  end  point.  Stated 
more  practically,  the  canonical  form  of  the  one-dimensional  equation 
for  elevated  receivers  is  inaccurate  unless  the  receiver  height  ex¬ 
ceeds  about  one-sixth  of  a  wavelength.  Note  that  this  problem  does 
not  occur  for  ground-based  receivers,  where  B^/Bn  =  0. 

To  better  quantify  the  end-point  error  caused  by  the  stationary- 
phase  approximation,  we  again  examine  Eqs.  (A. 6)  and  (A. 7).  The 
integration  just  outlined  is  valid  except  for  a  small  region  of  dimen¬ 
sion  ~z.  near  the  end  point.  Even  within  that  region,  only  the  term 
U  2 

proportional  to  ZQ/R2  in  Eci‘  (A. 7)  varies  so  rapidly  as  to  violate 
the  stationary-phase  approximation.  Thus,  we  can  apply  the  approxi¬ 
mation  over  the  entire  path  save  the  end-point  region,  where  the 
2 

Zq/r 2  term — which  stems  from  the  (l/k.R2) (BR^/Bn)  term  in  Eq.  (3) — 
requires  special  attention.  The  equation  that  results  is  identical 
to  Eq.  (A. 1) ,  except  that  the  "end-point  correction"  must  be  added  to 
the  right  side: 


end-point  correction 


(A. 17) 


where  the  upper  limit,  a,  should  be  taken  somewhat  larger  than  z^, 
but  much  smaller  than  x^;  and  9,  and  z  are  cylindrical  coordi¬ 
nates  centered  with  origin  on  the  ground  beneath  the  receiver.  For 

long  paths,  the  following  approximations  are  valid  within  the  end- 

1  /2 

point  region,  provided  zfl  <<  (Ax^)  : 
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«  R„  x. 


W(r2,  0)  «  W(xQ,  0,  0)  . 


With  these  approximations,  we  have 


(A. 18) 

(A. 19) 


a 

end-point  correction  W(Xg,  / 

•'0 


dr„ 


Z0r2 


ikR, 


(A. 20) 


The  integral  in  Eq .  (A. 20)  cannot  be  evaluated  in  closed  form 
because  of  its  finite  upper  integration  limit.  The  rapid  falloff 
of  the  integrand  for  r2  >  z^  permits  us  to  set  a  =  00  with  little  loss 
of  accuracy.  Then,  by  making  the  change  of  variable 


4  = 


V1  +  r2/z0 


we  f ind 


i 


t 


end-point  correction  «  W(x 


,°>/ 


«  ikzQC 
e 


di  =  W(xQ,  0)  E2(-ikzQ)  , 

(A. 21) 


where  E2  is  the  exponential  integral  whose  properties  are  given  by 
Abromowitz  and  Stegun  [1972],  By  using  the  power  and  asymptotic 
series  for  E2,  we  find 


W(xQ,  0)  , 


end-point  correction 


ikz. 


W(xQ,  0)  e 
Tkz~~ 


if  kzQ  «  1  ,  (A. 22) 


if  kzQ  »  1  ,  (A. 23) 


the  magnitudes  of  which  correspond  to  Eqs.  (18)  and  (19)  of  the  main 
text.  Note  that  adding  Eq.  (A. 22)  to  the  right  side  of  the  standard 
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one-dimensional  equation  for  elevated  receivers  [Eq.  (A.l)  or  Eq.  (5)] 
recaptures  the  continuity  of  the  fields  as  -*•  0.  Equations  (22)  and 
(23)  explain  why,  for  long  paths,  the  one-  and  two-dimensional  numer¬ 
ical  solutions  given  in  Sec.  Ill  agree  well  at  elevations  above  a  few 
kilometers  but  differ  substantially  near  the  ground. 

ELEVATION-ANGLE  ERROR  IN  ONE-DIMENSIONAL  APPROXIMATION 

In  addition  to  the  end-point  errors,  which  depend  on  the  ratio 
of  receiver  altitude  to  wavelength,  the  one-dimensional  equation 
suffers  from  an  elevation-angle  error  that  depends  on  the  ratio  z^/x^ 
of  receiver  altitude  to  pathlength.  It  is  intuitively  evident  that 
several  aspects  of  the  stationary-phase  integration  over  the  trans¬ 
verse  coordinate,  which  reduces  Eq .  (9)  to  Eq.  (A.l),  breaks  down  at 
large  elevation  angles.  Hufford's  [1952]  transformation  to  elliptical 
coordinates  fails  because  the  intersection  of  the  Fresnel  zone  ellip¬ 
soids  with  the  ground  strongly  deviates  from  an  ellipse  if  zq/xq  is 
Loo  large.  Similarly,  Field  and  Allen's  [1978]  treatment  encounters 
errors  because  as  zq/xq  increases,  the  stationary-phase  integration 
paths — the  intersection  of  hyperboloids  with  the  earth's  surface — 
are  less  accurately  approximated  by  the  lines  x  =  constant  (refer 
to  Fig.  A.l). 

One  way  to  quantify  the  elevation-angle  error  is  by  noting  that 
the  approx imat ion 


(A. 24) 


used  to  reduce  Eq.  (A. 6)  to  Eq.  (A.l),  fails  unless  zq/xq  is  small. 
When  the  receiver  is  on  the  ground,  the  following  relations  apply 
along  the  x  axis — which  closely  approximates  the  line  of  stationary 
points  for  long  propagation  paths: 

Godzinski  [1961]  shows  that  if  zq  >  *^xo»  higher  order  Fresnel 
zones  intersect  the  ground  in  a  manner  that  deemphasizes  the  contri¬ 
bution  of  points  near  the  receiver. 


The  exponents  in  Eqs.  (A. 6)  and  (A.1)  are  thus  invariant  between 
x  =  0  and  x  =  x^,  but  oscillate  rapidly  outside  that  range.  However, 
if  Zq  is  not  zero,  we  have  the  following  relations: 


ikz^ 

mrl  *  h  -  «„)  .  , 


ik(r^  +  R2  -  Rq)  ikzQ  , 


if  x0  -  X  >  z0  ; 


if  x0  -  x  «  z0  . 


(A. 26) 


If  kz0  <<  1,  no  substantial  oscillation  arises  in  the  region  0  s  x  S  xQ, 
and  the  upper  integration  limit  may  reasonably  be  taken  as  xR.  Recall 
however  that  the  end-point  error  becomes  substantial  if  kzQ  «  1. 

If  kz0  1»  Eq.  (A.  26)  shows  that  exp  [ikl^  +  R2  -  Rq)  ]  does 
oscillate  in  the  region  defined  by  xQ  -  <  x  <  xQ.  The  implication 

is  that  the  integrand  does  not  change  abruptly  from  nonoscillatory  to 
oscillatory  at  x  =  x^  as  it  does  for  a  nonelevated  receiver;  instead, 
change  occurs  over  a  transition  region  of  width  about  equal  to  zQ 
around  the  subreceiver  point.  Therefore,  unless 

z0/xQ  «  1  ,  (A. 27) 

the  assumption  that  the  integrand  slowly  oscillates  between  0  and  xQ 
will  be  invalid  over  a  significant  fraction  of  the  path. 

Although  the  criterion  given  as  Eq.  (A. 27)  derives  from  Field 
and  Allen's  [1978]  reduction  to  a  one-dimensional  equation,  it  applies 
to  any  stationary-phase  integration  over  the  transverse  coordinates 
[Hufford,  1952;  Johler  and  Berry,  1977],  That  point  is  borne  out  in 
Sec.  Ill,  where  the  disagreement  between  the  one-  and  two-dimensional 
numerical  results  is  shown  to  be  nearly  directly  proportional  to  the 
elevation  angle  z^/xq. 
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Appendix  B 

INTEGRAL  EQUATIONS  FOR  ISOLATED  IRREGULARITIES 
ON  A  UNIFORM  PLANE 


This  appendix  derives  Eqs.  (24)  and  (26),  which  are  more  con¬ 
venient  numerically  than  Eq.  (2)  for  small  terrain  irregularities. 

We  use  the  Green’s  function 

G(r2)  =  WH(r2WQ(r2)  (B.l) 

and  an  auxiliary  attenuation  function  W^ ,  defined  by 

<J»(P)  =  2W1(p)WH(p)^0(r)  ,  (B.2) 

instead  of  G  *  ^  and  ip  =  2Wi£q,  as  was  done  [Hufford,  1952]  to  derive 

Eq.  (2).  As  before,  <p  and  tp^  are  the  vertical  Hertz  potential  and  the 

Hertz  potential  in  free  space,  respectively.  W  is  the  attenuation 

n 

function  for  a  uniform  plane  given  by  Eq.  (22).  W^  therefore  accounts 

solely  for  impedance  or  topographic  irregularities.  As  shown  below, 

use  of  the  auxiliary  function  W,  instead  of  W  itself  (recall  W  E  W  W  ) 

1  1  H 

offers  computational  advantages. 

ISOLATED  LAKE  OR  ISLAND 

We  consider  a  plane  containing  no  topographic  features  so  that 

3r  /3n  is  everywhere  equal  to  zero.  The  refractive  index  is  n  (x,  y) , 
2  8 
which  equals  a  constant,  n^,  except  over  a  limited  region.  The  im¬ 
pedance  boundary  conditions  become 


_  JM 

3n  ’ 

g 


(B.3) 
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3G 

3n 


-  <P 


% 

0  3n 


ikG 

n 


(B.4) 


where  we  have  used  the  fact  that  3W„/3n  =  -ikW„/n  . 

H  H  g 

We  apply  Green's  theorem  to  obtain 


^  -  G(p) 


Zi  jji%  dy  [C<r2) 


3iKQ> 


3n 


3G( 

HQ) 


5Vl 

3n  J 


(B.  5) 


and  use  Eqs.  (B.3)  and  (B.4)  to  cast  Eq.  (B.5)  into  the  following 
form 


HKp)  =  2G(p)  + 


// 


dx  dy  G(r2)^(r1)  A(x,  y) 


(B.b) 


where 


A(x,  y)  =  ~  ~  ~~  (B.  7) 

n  n 
ng  g 

is  the  impedance  contrast.  Then,  using  the  fact  that 

ikr 

*0<r)  *  (B.8) 

and  inserting  Eqs.  (B.l)  and  (B.2)  into  Eq.  (B.6),  we  find 

ikro  lkro  ..  rr 

2Vp)Vr0>  ^  *  2VV  —  +  TjJ  d*  dy 

ik(r  +r  ) 

*  «1(Q)WH(r1)WH(r2)  A  e  .  (B.9) 


j 
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By  dividing  both  sides  of  Eq.  (B.9)  by  2W^(rQ)  exp  (itcr^/r^,  we 
obtain  Rq.  (24)  in  the  text.  Equation  (24)  is  unchanged  if  more  than 
uni-  hill  oi  lake  exists,  but  the  Integration  on  the  right  must  be 
performed  over  each  irregularity. 

ISOLATED  HILL 

The  derivation  of  the  integral  equation  for  in  the  presence 
of  an  isolated  hill  (or  hills)  proceeds  exactly  as  outlined  above  for 
an  isolated  lake  or  island.  We  note  that  here  the  refractive  index 
is  equal  to  n  over  the  whole  plane,  but  3^/30  does  not  vanish  every¬ 
where,  The  impedance  boundary  conditions  are: 


the  text. 
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Appendix  C 

ALGORITHM  FOR  NUMERICAL  SOLUTION  OF  TWO-DIMENSIONAL 
INTEGRAL  EQUATION 


This  appendix  gives  the  algorithm  used  to  obtain  numerical  solu¬ 
tions  of  Eqs.  (24)  and  (26),  which  can  be  written  in  the  form 


W(x,  y) 


dy' 


k(x,  y. 


y’)W(x\  y')  . 


(C.l) 


If  terrain  features  are  localized,  the  kernel  of  Eq.  (C.l)  can  be 
assumed  to  vanish  outside  some  bounded  region.  The  solution  consists 
of  first  solving  Eq .  (C.l)  to  find  W  within  this  irregular  region, 
and  then  integrating  over  the  region  to  obtain  W  on  the  uniform  por¬ 
tion  of  the  plane. 

SOLUTION  WITHIN  REGION  OF  IRREGULARITY 

We  divide  the  isolated  region  into  a  square  grid  as  shown  in 

Fig.  C.l,  where  n  and  n  are  the  number  of  squares  in  the  x  and  y 
b  x  v 

directions.  Equation  (C.l)  then  becomes 


n  -1  n  -1 
x  y 


W(x,  y)  =  1  + 


£  £  //-  *• 


K(x,  y,  x’,  y’)W(x’,  y ’ )  ,  (C.2) 


i-0  j=0 


Sij 


where  on  the  square  s^,  centered  at  x^,  y j , 


As  ^  ^  .As 

xi  -  -y  ;£  x’  £  xi  +  y  , 


As  ,  .As 

71  /  '  v  C  v)  ’  • 


(C.3) 


Fig.  C.l--Parti 
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Here  As  is  the  length  of  the  side  of  a  square  taken  to  be  small  enough 
that 


W(x,  y)  «  W(xi>  y  ^ ) 


(C.4) 


on  each  square.  We  can  then  rewrite  F.q.  (C.2)  as 


n  -1  n  -1 
x  y 


W(V  yt)s!l+ES  W(V  V  ffdK 

1=0  j=0  JJ 

s .  . 

ij 


dy'  K(xk,  y j ,  x’ ,  v') 


(C.5) 


By  defining 


\uj  ffdx'  dy'  K(V  >V  x’*  y’: 


’ij 


the  solution  to  Eq.  (C.l)  can  be  expressed  as 


(C.<>) 


n  -1  n  -1 
x  y 


+  £  X)  W<V  V3 


)M, 


kfcij 


W(xk,  yj  = 


i=0  j=0 

i^k 


1  -  M. 


(C.  7) 


kfck£ 


The  integration  in  Eq.  (C.6)  is  performed  with  standard  quadrature 

methods.  We  use  a  higher  density  of  quadrature  points  when  (i,  j)  = 

(k,  H)  to  accommodate  the  proximity  of  the  singularity  in  K  at  r^  =  0. 

Equation  (C.7)  is  a  system  of  n  *  n  linear  equations,  which 

x  y 

we  solve  iteratively  using  the  Gauss-Siedel  method  with  each  W  initially 

set  equal  to  one.  To  evenly  distribute  errors  throughout  the  grid, 

we  begin  each  iteration  at  a  different  point.  Arbitrary  convergence 

criteria  can  be  used.  Computation  time  varies  as  (n  x  n  )^,  so  care 

x  y 

must  be  exercised  in  selecting  the  parameter  As.  In  most  problems, 
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che  optimum  value  for  As  depends  on  wavelength,  rather  than  on  terrain 
features.  Computation  time  then  varies  as  the  square  of  the  area  of 
the  perturbed  region,  and  limits  the  size  of  irregularities  that  can 
be  analyzed  at  reasonable  cost.  For  a  frequency  of  100  kHz,  we  found 
the  practical  limit  on  size  imposed  by  computational  cost  to  be  about 
10  km  x  10  km. 

SOLUTION  OUTSIDE  REGION  OF  IRREGULARITY 

Once  W  has  been  found  inside  the  irregular  region,  its  value  out¬ 
side  is  calculated  from  Eq.  (C.5),  using  weli-known  quadrature  methods. 
Here  W  is  given  by  a  double  integral  over  this  region,  rather  than  by 
a  solution  to  a  two-dimensional  integral  equation.  Also,  the  singu¬ 
larity  at  =  0  causes  no  difficu  lty. 
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